N94- 


A/ /-S 

27287 


Preprint, AIAA Atmospheric Flight Mechanics Conference. Hilton Head, SC, USA, August 10-12 1992 


2e 757/ 


A Data Fusion Algorithm for Multi-Sensor Microburst Hazard Assessment 


Craig R. Wankc* and R. John Hansman** 
Department of Aeronautics and Astronautics 
Massachusetts Institute of Technology 
Cambridge, MA USA 



> 


Abstract 

A recursive model-based data fusion algorithm for 
multi-sensor microburst hazard assessment is described. 

An analytical microburst model is used to approximate the 
actual windfield, and a set of ‘"best” model parameters are 
estimated from measured winds. The winds corresponding 
to the best parameter set can then be used to compute 
alerting factors such as microburst position, extent, and 
intensity. The estimation algorithm is based on an 
iterated extended Kalman filler which uses the microburst 
model parameters as state variables. Microburst state 
dynamic and process noise parameters are chosen based on 
measured microbursl statistics. The estimation method is 
applied to data from a time-varying computational 
simulation of a historical microburst event to demonstrate 
its capabilities and limitations. Selection of filter 
parameters and initial conditions is discussed. 

I Computational requirements and datalink bandwidth 
i considerations are also addressed. 

InlrQductifln 

Low altitude wind shear has been a major cause of 
fatal aviation accidents in the U.S. 1 The localized intense 
downdrafts known as microbursts are the most dangerous 
form of wind shear, and pose a serious hazard to aircraft 
during takeoff or approach. In a typical microburst 
encounter, an aircraft first encounters a performance- 
increasing headwind. This is followed by a downdraft and 
a rapid transition from headwind to tailwind, which 
produce sharp losses in altitude and/or airspeed. 

Several systems for detection and measurement of 
microburst hazards are currently nearing the operational 
stage. Effective ground-based systems such as Terminal 
Doppler Weather Radar (TDWR) and the extended Low 
Level Wind Shear Alert System (LLWAS) are entering the 
deployment phase. TDWRs will be located at 47 major 
airports, and detect microbursts primarily by measuring 
the surface wind velocity component radial to the radar and 
identifying areas of radial shear. 2 - 3 LLWAS is a network 
of anemometers which measure horizontal windspeed and 
direction around the airport surface, and detect wind shear 
events from differences in wind speed and direction 
between sensors. 4 Airborne reactive wind shear alerting 
systems, currently in use, detect microburst penetration by 
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comparing inertial and air data system measurements to 
compute the local winds. Several types of airborne 
forward-look sensor technologies are also under 
development, including infrared radiometry, Doppler radars 
and Doppler lidars. 5 Infrared systems measure the drop in 
temperature associated with the air in the center of a 
microburst, while Doppler radars and lidars measure wind 
velocities along the flight path ahead of the aircraft. In 
addition to new sensor developments, the development of 
digital air- ground datalink capabilities such as the Mode-S 
beacon system will allow microburst alert information to 
be exchanged between air and ground-based systems 
(Figure l). 6 

As new detection systems become operational, it will 
become likely that more than one sensor system will be 
available in a given situation. Also, each of the 
aforementioned sensor systems has some geometrical 
observability problems. For example, both ground-based 
and airborne Doppler radars and lidars can only measure 
wind velocities radial to the sensor, not vertical winds. 
The aviation hazard posed by a microburst, however, is 
due to both horizontal wind shear and downdrafts in the 
microburst core. Therefore, a technique for combining 
data from different systems with different measurement 
characteristics could improve estimates of microburst 
hazards and aid alert generation. 

The goal of this “data fusion” process is to provide a 
microburst detection and hazard assessment capability 
which is significantly better than that which can be 
achieved using a single sensor. The data fusion algorithm 
must provide appropriate information for alert generation, 
in a timely fashion, and be feasible with regard to the 
available air-ground datalink bandwidth and computational 
capabilities. Previous work at MIT has focused on 



Figure 1: Advanced Microburst Detection 

and Alerting Systems 
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definition of appropriate information for microburst alerts. 
This work has included analytical studies to determine 
appropriate microburst hazard criteria, 7 and piloted part- 
task simulator studies to determine the information 
requirements of the flight crew when faced with a 
microburst alert situation. 8 9 The conclusions of these 
studies were that estimates of microburst location, extent, 
and approximate intensity are required for alerting 
purposes, and that microburst intensity can be quantified 
well using criteria which relate to the expected aircraft 
energy loss due to the microburst windfield. Therefore, a 
good data fusion algorithm should be able to compute the 
position and extent of a microburst, and to estimate 
intensity including contributions from both horizontal 
wind shear and downdrafts. 

The data-fusion process can be done on a number of 
levels. One approach is “to merge the final products of the 
sensor systems to produce improved alerts. For example, 
product-level techniques have been used to integrate 
TDWR and LLWAS information 10 and to determine the 
probability of hazardous wind shear given a wide range of 
evidence.! 1 Another approach is to integrate sensors on 
the data level. The data-level approach is more complex, 
due to the large volume of data produced by several wind 
shear sensor systems. However, if correctly implemented, 
observability problems due to poor sensor geometry can 
be alleviated. Data from multiple sensors can be 
combined to form a "super sensor 1 ' which has improved 
sensing geometry. The technique proposed in this paper 
is a model-based data-level approach which attempts to 
gain this observability benefit without prohibitively large 
computational or data transmission requirements. 


Why Use a Model? 

Representation of the actual microburst windfield 
with an analytical model has two major advantages. 
Firstly, once the model has been "fitted" to the windfield, 
the windfield can be approximated by the values of the 
model parameters. Thus, if the model represents the 
windfield well enough, the measured information (a large 
data set) can be encapsulated in a small set of "best-fit” 
model parameters. Since it is impractical (or at least 
undesirable) to transmit raw data between aircraft and 
ground-based systems, this is an important advantage. 

Secondly, an analytical model can include additional 
information which can be used to infer quantities which 
cannot be directly measured, such as inferring vertical 
velocities from radar-measured radial velocities. 

Analytical models can be designed to satisfy basic fluid 
dynamic relationships such as mass continuity, and can be 
adjusted to reflect results of microburst field 
measurements. 

Analytical Microburst Mode! 

The analytical microburst model used in this work 
was developed at NASA Langley Research Center initially 


by Oseguera and Bowles, 12 and later improved by 
Vicroy. 1 3 The Osegucra-Bowles-Vicroy (OBV) model 
uses shaping functions to generate an axisymmetric 
flowficld which satisfies the mass continuity equation and 
is representative of the major characteristics of measured 
microbursts. Sample winds for a constant -altitude path 
through the model windfield are shown in Figure 2. The 
horizontal winds exhibit the classic microburst 
characteristic of a headwind increase, followed by rapid 
shearing to a tailwind. The vertical wind plot shows a 
downdraft in the microburst center and smaller updrafts at 
the edges. 

The microburst w inds are uniquely defined by a set of 
five parameters and three empirically -adjusted constants. 
For this stu dy, a simple ambi ent wind (4 additio nal 
parameters) was added to the microburst windfield. The 
model parameters are summarized in Table 1. The total 
winds are given by non-linear, smooth, differentiable 
functions of the parameters and a given (x,y,h) position as 
follows: 

U ~ LImicro burst + Uo + Uhh (1) 

V = Vmicroburet + Vq + Vhh (2) 

W = W m j cro bursl (3) 

U, V, and W are the Eastward, Northward, and vertical 
wind velocities; h is the altitude above ground level. 
"Microburst" quantities are functions of position and of 
the first five parameters in Table 1; these functions are 
summarized in Appendix A. 

The OBV model is axisymmetric, but naturally 
occurring microbursts are often asymmetric. 14 In 
addition, multiple microbursts have been observed to 
occur close together and interact. To handle these cases, 
the model was extended to allow multiple interacting 
microbursls. For each microburst, another set of five 
microburst model parameters (the first five in Table 1) can 
be added. It is assumed that the ambient wind will be 



x (m), distance along flight path 


Figure 2: Microburst Model Windfield. Sample 

winds for a constant-altitude flight path through the center of a 
simulated microburst. 
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Table 1: Modified Oseguera-Bowles-Vicroy 

Microbursl Model Parameters 


Parameter 

Description 

*0 

X-coordinate (Hast) of microburst center (m) 

yo 

Y-coordinate (North) of microburst center (m) 

U m 

Maximum horizontal outflow speed (m/s) 

R d 

Radius of maximum outflow (meters) 

Zn, 

Altitude AGL of maximum outflow (meters) 

u 0 

Eastward ambient wind constant component (m/s) 

u h 

Eastward ambient wind altitude gradient (m/s/m) 

v 0 

Northward ambient wind const, component (m/s) 

v h 

Northward ambient wind altitude gradient (m/s/m) 


roughly constant throughout the x-y space of interest, i.e. 
near the airport, and so only one set of ambient wind 
parameters is used. The winds from each model 
microburst are summed to get the overall model windfield; 
this superposition does not violate mass continuity. In 
the simulation study below, when a “two-microburst” 
model is referred to it does not necessarily indicate that 
there are two microbursts being detected, but that two 
superimposed model microbursts are being used to 
simulate a complex microburst windfield with more than 
one area of high downdraft speed. 

Model-Based Multi-Sensor Data Fusion 

Given a suitable model, the fusion problem reduces to 
estimation of the "best" set of model parameters based on 
all available wind measurements. Once the be st 
parameters have be en estimate d, alerting. f actorslmteDiil^ 
extent etc.) can be derived from the ana ly tic al mo del 
windfield corresponding to th e estimated model 
parameters. 

This estimation procedure must satisfy several 
constraints to be practical. The estimation algorithm 
must be recursive, to handle new measurements as they 
become available. It must also account for time variation 
in the model parameters, since microbursts are dynamic 
phenomena with short lifetimes on the order of 15 
minutes and sensor measurements will be taken at 
different times. It should also be probabilistic, to take 
advantage of microburst statistical characteristics from 
past field studies. A Kalman filter approach is proposed 
to satisfy these requirements. 

Iterated Extended K alman Filter Algorithm 

estimation Problem Structure 

Kalman filtering techniques require a state-space 
dynamic model of the system and a relationship between 
system parameters and measured quantities. In this case, 
we would like to estimate analytical model parameters 
which best describe the microburst from measurements of 
the winds. The analytical microburst model parameters 
are therefore used as the filter state variables x(t). It was 
assumed that the time evolution of the microburst 


parameters can be adequately modeled by a linear, time- 
invariant, continuous-time system: 

x(t) = A x(t) + B u(t) + L w(t) (4) 

Deterministic inputs to the system are represented by 
u(t), and w(t) is a white Gaussian process noise input. 

The A, B, and L matrices define the dynamic model; they 
will be discussed below. Since the state variables x(t) are 
the analytical model parameters, they are related to the 
wind measurements through the analytical model wind 
equations. The resulting non-linear discrete-time 
measurement equation is: 

Zk = hk(x(tk)) + vk (5) 

where the measurement equations hk are simply the 
wind equations from the analytical model, and Vk 
represents measurement noise. The state vector, x, and 
error covariance matrix, P, for a single downdraft case are 
defined as follows: 

x=[x 0 yo U m Rp z™ Uo U h V 0 V h ] T (6) 

P = e[(x - x)(x - xf] (7) 

where x is the current parameter estimate. The 
microburst eastward core location, xo, should not be 
confused with x, the slate vector. Process noise, w(t), and 
measurement noise, v k , are white and gaussian with the 
following characteristics: 

E [(L w(t)) (L w(x)) T ] = L Q(t) L t 8(t - z) (8) 

E [vk Vk] = R k (9) 

The aim of the filter is to produce the state estimate x 
which minimizes the error covariance P. Since the 
measurement equation is non-linear, this cannot be done 
with a standard Kalman filter algorithm. An extended 
Kalman filter (EKF) approach was chosen. The structure 
and principal equations for the EKF are briefly described 
below, based on the formulation given in Ref 15. The 
filtering algorithm for discrete-time measurements is a 
two step process: (1) apply the system dynamic model to 
propagate the state estimate and state estimation error 
covariance between measurements, and (2) update the 
estimate when new measurements arrive. 

Estimate Propagation: Microburst Dynamic 

Modeling 

For linear, time-invariant, continuous-time system 
dynamics the propagation of the state estimate and 
estimation error covariance between measurements is 
governed by: 

x(t) = A x(l) + B u(t) (10) 




P(t) = A P(t) + P(t) A T + L Q(t) L t 


( l *) Kkj = Pk Hk (xk.i)[Hk(xkj) Pk Hk(xk,i) + R k ] (15) 


The A, B, L, and Q matrices define the microburst 
time-evolution dynamics. Since (he analytical model is 
lime -invariant, these parameters must come from another 
source. Unfortunately, there is no simple time-varying 
analytical model available. However, measured 
microburst statistics can be used to approximate some 
dynamics" For example, microburst radial extent tends to 
increase steadily throughout the microburst lifecycle. 
Analysis of data from Colorado microburst 
measurements 16 - 17 indicates that the change in radial 
extent vs. time can be approximated by a constant bias (a) 
with additive white zero-mean gaussian noise (n): 


R p = a + n(t) 


( 12 ) 


where a = 0. 102 km/m in and the noise term has a 
standard deviation of 0. 15 km/min. The constant bias is 
treated as a deterministic input, and the noise term leads to 
a value for one element in the Q matrix. Similar 
modeling may be possible for some of the other state 
variables. For example, motion of the microburst core 
(Xo,y D ) may be related to the ambient wind parameters, 
which would lead to non-zero entries in the A matrix. In 
the simulations discussed below, the A matrix was 
assumed to contain all zeros. The B and Q matrix 
elements were set based on statistical information where 
possible, and from engineering judgement when no 
statistical information was available. Further research on 
microburst dynamics is currently in progress. 

Since the time behavior of the microburst parameters 
is not well modeled, significant process noise is required. 
The use of process noise to compensate for modeling 
deficiencies is similar to the well-known technique of 
applying a “forgetting factor” to older data in a batch least- 
squares formulation. In any case, these simple dynamics 
lead to sparse A, B, L, and Q matrices, and the 
propagation step in the filter requires little computation. 

Incorporating Measurements 

When new measurements are taken, the estimate is 
updated. The non-linear measurement equation, however, 
makes the update process difficult. The formulation 
presented here is based on the extended Kalman filter 
update with the addition of a local iteration procedure to 
reduce the effects of the measurement non-linearities, 15 - 18 
At time a local iteration (over i) is performed. The i* 

parameter estimate at time t*, xj^, is updated with the 
following expression: 

Xk,i+i = Xk + Kkj [zk - hk(x£,J - Hk(x^) (xk - XkJ] (13) 
xic.o^xk, i = 0,l,... (14) 

The Kalman gain, K, is ordinarily computed from: 


and H k is the locally linearized measurement matrix: 


Hk(x) = 


ahk(x) 


dx 


(16) 


In the above expres sions, Xk and Pk indicate the 
propagated estimate and error covariance at time tk (prior 

to updating), while Xk and Pk indicate the updated 
estimate and covariance based on the measurement The 
local iteration is repeated until the scaled norm of the 
parameter estimate does not change significantly. After 
the new estimate has been produced, the updated error 
covariance matrix is computed using values from the final 
iteration step: 

Pk = [i - Kkj H k (xu)] Pii (17) 

Some simple testing, in which winds generated 
directly from the OBV model were “identified” using this 
algorithm, indicated that the iterated filter results in 
significantly better estimates than the standard EKF; this 
has also been found by other investigators. 19 A 
probabilistic interpretation of this iteration based on 
Bayesian maximum likelihood estimation is given in Ref. 
18. 


One difficulty with the above updating algorithm is 
that there may be large numbers of measurements 
available at a single time step (as in TDWR data, for 
example), and the computation of the Kalman gain (Eqn. 
15) requires inversion of an r-by-r matrix, where r is the 
number of measurements. The number of computations 
required to do this scales as r 3 . In a linear filter, a large 
batch of measurements can be treated as a series of 
sequential scalar measurements (occurring at infinitesmal 
time spacing) without loss of information, thereby 
avoiding this problem. When the measurement equation 
is non-linear, the measurements cannot be incorporated 
sequentially without losing a significant amount of 
information. Therefore, an alternate form of the gain 
computation is required. When the number of 
measurements exceeds twice the number of states, and the 
measurement noises arc independent (diagonal Rk) it is 
more efficient to use the ” information form" of the gain 
computation: 

(Pk'JMPk'J' + HjRi'Hk (18) 

MPkTtfRk 1 09) 

This form can be readily applied to the iterated EKF 
update described above. Although the covariance update 
must now be done inside the loop, the required matrix 
inversion is only n-by-n, where n is the number of states 
(model parameters). The computational requirement now 
scales linearly with r and cubicaily with n. In the 
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simulation cases below, where r > 100 and n = 9 or 14, 
this form was found to be much more efficient. 

Multiple-Microburst Form 

As discussed above, several model microbursts can be 
superimposed to simulate a more complex windfield. In 
this case, 5 new states are added for each additional model 
microburst. For i microbursls, the full stale vector, x, is 
defined as: 

x-tiu W--.xl.iUo Uh Vo Vhf (20) 

where: 

Xmb.i^fxo.i yo.i Um,i Rp.i Z mt i] (21) 

In the simulations discussed below, one- and two- 
microburst forms are used. 

Initialization 


point). This polygon encloses the entire performance- 
decreasing portion of the microburst. The ability of the 
model-based algorithm to define this hazardous region can 
then be evaluated by comparing the extent polygon A of 
the truth windfield to the extent polygon B of the 
analytical model windfield corresponding to the estimated 
parameters: 

= ( 22 ) 
AuB 

This quantity has a maximum value of 1 (for an exact 
match) and falls off for both underestimation and 
overestimation of the extent boundaries (Figure 4). Core 
position errors are also reflected, since the model extent 
polygon is then laterally translated with respect to the true 
extent polygon. This quantity is a function of altitude, 
but the dependence was found to be very weak and only 
results for a single altitude are presented in this paper. 


This algorithm can incorporate multi-sensor data, 
given that the microburst has previously been detected. 

The assumption is made that a single sensor has detected 
the event and has produced an initial parameter estimate 
and associated error covariance. The initialization 
algorithm therefore depends on the measurement 
characteristics of the initial sensor. The general process, 
however, is the same for all sensors. Quantities that can 
be directly measured are estimated from the initial data set, 
and quantities which are unobservable are initialized using 
statistics derived from microburst field studies. 

For example, if a TDWR initially detects a 
microburst, estimates of maximum outflow speed (U m ), 
outflow radius (R p ), and core position (x 0 ,yo) can be 
derived from the radial flowfield measured by TDWR. The 
outflow depth is an unobservable parameter, and must 
be initialized from statistics. Outflow depth statistics 
have been measured for 26 Colorado microbursts, 20 and 
the mean altitude of maximum outflow velocity was 
found to be 109 meters. This value was used to initialize 
the Filter for the simulations discussed below, in which 
TDWR was always assumed to make the initial detection. 
The initial covariance matrix was diagonal, and values 
were chosen based on sensor resolution criteria or 
statistics where possible. 

Simulation Results 

Figures of Merit 

As mentioned above, the important quantities for 
alerting purposes are position, extent, and approximate 
intensity. The “effectiveness” of the proposed algorithm 
can be defined in terms of its capability to produce these 
quantities. Therefore, two figures of merit were defined. 
The First concerns position and extent. Given a center 
point, an “extent polygon” can be drawn for a microburst 
windfield (example shown in Figure 3). 14 The vertices of 
the polygon correspond to the points of maximum radial 
outflow speed (measured radially outward from a center 



Figure 3: Microburst extent polygon. Horizontal 

windfield with outflow extent polygon superimposed. 



Figure 4: Extent Figure-of-Merit. Pictorial 

representation of Eqn. 22. The cross-hatched area is A o B 
and the sum of the cross-hatched and striped areas is A B - 
Note that the model “polygon” B is a circle for the single- 
microburst axisymmetric OBV model. 
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Microbum intensity was defined in terms of "F- 
factor", proposed by researchers at NASA Langley 
Research Center, which is based on the impact of a 
microburst windfleid on the total energy (lcinetic plus 
potential) of the aircraft. 5 It is a measure of the loss of 
potential rateof-climb (or loss of effective thrust-to- 
weight ratio) due to the immediate windfield. It is 
dependent on the time rate of change in the aircraft frame 
of the tailwind velocity, the vertical wind velocity, and the 
aircraft airspeed. Positive values of F indicate a 
performance-decreasing situation, and negative values 
indicate a performance -increasing situation. As typical 
transport-category aircraft in landing or takeoff 
configuration have excess thrust-to-weight ratios between 
0.1 and 0.15, an encounter with an F-factor in excess of 
that value would compel the aircraft to descend and is 
therefore hazardous. 

g V 

F, a point measurement, needs to be averaged over a 
distance to give a useful indication of aircraft hazard. Past 
work has determined that F averaged over 1 km of the 
aircraft flight path yields a good hazard estimate. 7 For 
evaluation purposes, however, it is desirable to assign a 
single hazard number to a microburst rather than one for 
each possible flight path through it. Therefore, for this 
work, the hazard number was defined as follows: (1) 
compute 1 km average F-factors for a large number of 
parallel constant-al t it ude paths through the microbunt. (2) 
average the resulting values 500 m laterally across flight 
paths, and (3) pick the largest averaged F-factor as the 
hazard value. 

This value depends strongly on the direction of the 
flight paths along which F is evaluated. In the simulation 
results presented below, averaged F-factors will be 


presented for either eastbound or northbound flight paths. 
In addition, F depends on altitude, and results will 
therefore be presented for several altitudes. For alerting 
purposes, however, it would be necessary to assign a 
single intensity value to a detected microburst, for : 
example the largest value (over all directions) below a 
specified maximum altitude. 

Simulated Microburst True Windfield 

The windfield data used to evaluate the estimation 
algorithm was generated by the Terminal Area Simulation 
System (TASS). 21 It is a highly detailed computational 
simulation of a complex multiple microburst event which 
occurred at Denver-Stapleion airport on July 1 1. 1988. 
This event caused one near accident and a total of five 
aircraft to make missed approaches. 22 Windfield data 
from five times during this event was available, with a 
horizontal spacing of 200m and a vertical spacing of 
approximately 80m. For the following analyses the 
largest microburst in the event was selected (Figure 5). 
The horizontal windfield has a classic microburst outflow 
pattern. However, the vertical wind contours show some 
complex structure as indicated by two separate regions of 
high downdraft, neither of which correspond to the 
apparent horizontal windfield center (marked with an X). 
This rather complex event was chosen to test the 
estimation algorithm in a challenging but realistic 
situation. 

OBV Model Best Fit 

The first step for algorithm evaluation was to 
determine the ability of the analytical model to match the 
important characteristics of a microburst windfleid, 
namely the Figures of merit defined above. This was done 
with a deterministic non-linear batch least-squares 
optimization algorithm, similar to that used in Ref. 23 to 
model microburst winds with vortex rings. The “truth" 



Figure 5. TASS-simulated windfield for 11 July 1988 microburst event at DEN. At left is a vector plot of 
horizontal winds: maximum velocity shown is 18.7 m/s. At right: vertical windspeed contours. Altitude shown is 177m AGL. 
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winds were taken from one time-step of the TASS model 
(shown in Figure 5). The data volume included three-axis 
winds at three altitudes with 200m lateral spacing for a 
total of approximately 5000 data points. The 
optimization procedure was a constrained version of the 
standard Gauss-Newton method, 24 and found the model 
parameters which minimized the mean square wind error, 
J: 


N 

J = -L^ ej^p) ej(p) ; e(p) = V, mth - V mo( k<p) (24) 
N ]-i 

where N is the number of total data points and V is 
the vector of all wind points, including East, North, and 
vertical components at all (x,y,z) locations. This 
procedure was done with a single-microburst model (9 
parameters) and with a two- microburst model (14 
parameters). The resulting parameters for the single- 
microburst case are given in Table 2. Note that the 
approximate radius of this microburst is 1700 meters, and 
the maximum outflow speed is approximately 18 m/s. 
The ambient wind magnitude is small in this case. 

Table 2: Single-mlcroburst least-squares 

parameter fit results 


xo (m) 

9528 

U 0 (m/s) 

0.9 

yo (m) 

-5047 

Uh (m/s/m) 

-0.001 

U m (m/s) 

17.8 

V 0 (m/s) 

0.5 

R d (m) 

1717 

Vh (m/s/m) 

-0.002 

An ( m ) 

68.2 

1 


The single-microburst fit produced an extent figure of 
merit of 0.92. The two-microburst Fit result was slightly 
lower, at 0.85. As seen in Figure 5, this microburst was 
fairly axisymmetric in extent, so these good results are 
not surprising. However, plots of area-averaged F-factor 
looking Eastward and Northward for three altitudes (Figure 
5) reveal that the microburst is not symmetric in 
intensity. As indicated by the 'TASS windfield" points in 
Figure 6, the F-factors are larger when looking northward 
through the microburst than when looking eastward. This 
is due to the vertical wind distribution (Figure 5, 
righthand plot) which has multiple regions of high 
vertical windspeeds. For this reason, the single - 
microburst Fit produces a single broad region of somewhat 
weak vertical winds in an attempt to globally match the 
windfield, and the result is that intensity is underestimated 
in both directions. The two-microburst Fit, on the other 
hand, succeeds in matching the vertical windfield well and 
duplicates the intensity of the TASS windfield well in 
both directions. 

For alerting purposes, both model windfields 
adequately represent the actual extent; however, the single- 
microburst model underestimates the intensity somewhat. 
The results of previous work, however, indicate that 
highly accurate intensity estimates are not critical for alert 
generation 9 Based on these results, and similar results 
obtained using TASS windfields from another microburst 
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event, the model was judged to be acceptable for 
estimation purposes. 


Iterated EKF Update Simulation 


The next step in algorithm evaluation was to 
determine if the iterated EKF update procedure was capable 
of taking sensor data (as modeled by small subsets of the 
entire windfield) and producing reasonable extent and 
intensity estimates. The TASS simulated winds were 
again considered to be the ‘'truth'’ winds, and simulated 
sensor data subsets were taken from them. Assuming that 
the windfield was frozen in time (or alternatively, no time 
has lapsed between measurement sets), different 
combinations of sensor data were used sequentially to 
update the current estimate. Three sensors were considered 
in this way: (1) TDWR data, (2) winds measured from the 
aircraft, using inertial and air data measurements (referred 
to henceforth as INS data), and (3) airborne Doppler radar 
(ABDR) data. 

For TDWR and ABDR data, it was assumed that the 
sensor was far enough from the microburst that radial 
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wind measurements could be considered parallel to each 
other, and that the antenna tilt angle was horizontal so 
that all data was taken at the same altitude. For example, 
for an eastward -looking radar, the U-components of the 
TASS windfield at a single altitude became the working 
data set. TDWR measurements were taken at an altitude 
of 82m AGL (the lowest TASS data altitude) for both 
eastward- and northward-looking cases, and ABDR 
measurements were taken from 177m and 283m AGL 
TASS data. In all cases, gaussian zero-mean white noise 
with a standard deviation of 1 m/s was added to the ‘‘truth” 
data to simulate measurement noise based on TDWR 
accuracy specifications. 25 All radar data sets were taken at 
400m range and azimuth resolution; this is poorer than 
the resolution of operational radars, but reduced the 
computation time required to run the simulations. 

Aircraft winds (INS data sets) consisted of 3- 
component winds along a straight flight path at constant 
altitude. Four 200m resolution INS data sets were 
defined, including eastward and northward flight paths at 
177m and 283m AGL. AH paths passed through the 
approximate center of the windfield, as marked in Figure 
5. The measurement noise standard deviation used for 
aircraft wind measurements was 1.4 m/s. 

For simulation purposes, it was assumed that TDWR 
made the initial microburst detection. Therefore, the first 
step was to initialize the filter as previously described, and 
then apply the iterated EKF update to incorporate the 
TDWR measurement. The resultant parameter estimate 
and error covariance were saved. Then the estimate was 
updated by incorporating either an INS data set or an 
ABDR data set, starting with the saved parameter estimate 
and covariance matrix. Twelve total sensor fusion cases 
were tested with both one-microburst and two-microburst 
versions of the filter. 

Single-Microburst Filter 

For all cases tested, the iteration procedure used in the 
update converged in 3 to 5 iterations. Results for four 
representative cases are presented here; 

(1) Initialization only: Eastward-looking TDWR 
measurements alone (denoted TDWR-E) 

(2) The results of ( I) were updated using a sequence of 
eastbound aircraft-measured winds taken at an altitude of 
I77m AGL (denoted INS-E) 

(3) The results of (1) were updated using a sequence of 
northbound aircraft- measured winds taken at an altitude 
of 177m AGL (denoted INS-N) 

(4) The results of (1) were updated using northward- 
looking airborne Doppler radar data at 177m AGL 
(denoted ABDR-N) 

The extent results are again good (Table 3), and 
illustrate the effect of fusing data from sensors with 
different measurement geometries. The extent figure-of- 
merit for TDWR-East is 0.85, and does not improve when 


Table 3: Extent figures-of-merit: 1-microburst 

time-invariant data fusion. 


TDWR-E 

alone 

TDWR-E 
+ INS-E 

TDWR-E 
+ INS-N 

TDWR-E 
+ ABDR-N 

0.853 

0.853 

0.911 

0.917 


an eastward path of INS data is incorporated. However, 
when northbound INS data or northward-looking airborne 
radar data is incorporated, the extent figure-of- merit 
increases to the .9T to .92 range. Since the microburst is 
not exactly symmetric in extent (it is slightly larger in the 
north-south direction), incorporation of north ward- looking 
data increases the radius parameter in the OBV model to 
cover more area. This is equal to the performance 
achieved by the least-squares fit computation. 

The effect of multi-directional data is also visible in 
the intensity results (Figure 7). As with the least-squares 
results, it is clear that the single-microburst model cannot 
match intensity with the complex windfield of this 
microburst. The TDWR-alone result is low, and 
incorporating an eastbound path of INS data actually 
lowers the estimate; this is because the path does not 
cross both regions of high vertical windspeed. 

Incorporating a northbound path of INS data or the ABDR 
data improves the estimate significantly at the higher 
altitudes, from which the INS and ABDR data are taken. 

Two-Microburst Filter 

The two-microburst version of the filter involved 
significantly more computation, since in general more 
iterations were required than for the single-microburst 
filter. Also, some cases did not converge consistently and 
required adjustment of the initial parameters. However, 
when the two-microburst filler did converge, the results 
were good. Extent figures-of-merit were between 0.85 and 
0.90 for all cases. Figure 8 shows eastward inlensity 
values for the algorithm applied to three cases: 
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(1) northward-looking TDWR data alone (TDWR-N) 

(2) TDWR-N updated with eastward aircraft -measured 
winds taken at 177m AGL (INS-E) 


Case 1. A sequence of eastbound aircraft-measured 
winds is downlinked to the ground and incorporated 
along with a second set of TDWR data 

Case 2. An aircraft traveling northbound receives the 
previous TDWR estimate and updates using an airborne 
Doppler radar 

Case 3. An aircraft traveling eastbound receives the 
previous TDWR estimate and updates using an airborne 
Doppler radar 

At the third time step (+4 minutes), the parameter set 
is passed to the ground and an update is done using 
another set of TDWR data. The data sets were derived in 
the same way as for the lime- invariant cases, and the 
estimate and error covariance were propagated between 
measurements as described above. In all cases, the single- 
microburst form of the Filter was used. 

The extent figures-of-merit (Figure 9) are fairly good 
(> 0.82) through the first two times, but are slightly 
lower (0.76) in the third time step. This is due to the 
distorting effect of an adjacent, weaker microburst on the 
shape of the primary microburst. The axisymmetric 
model used in the filter has difficulty representing this 


(3) TDWR-N fused northward-looking airborne radar 
data at 283m AGL (ABDR-N). 

Although the TDWR is looking north and the 
intensity values shown are for eastbound paths, the results 
for TDWR data alone are fairly good. There is some 
overestimation at high altitudes. Inclusion of ABDR data 
with the same look angle as the TDWR (northward) 
improves the results slightly. As expected, inclusion of 
the eastbound INS data provides a second measurement 
direction and produces the best intensity estimates. 

Full Iterated EKF Simulation 

The third part of the algorithm evaluation was to 
include the microburst dynamic model (the propagation 
part of the niter) and apply the technique to time-varying 
data. For this analysis, data was taken from three different 
times in the evolution of the 7/11/88 microburst event. 
The three data sets were spaced two minutes apart, where 
the middle data set corresponds to the time-invariant data 
set used in the previous section and corresponds to the 
time at which the microburst was strongest. 

The time spacing for this data was larger than desired, 
since TDWR data is updated at 1 min intervals and 
airborne radar data would be available even more 
frequently. However, it was still possible to construct 
illustrative examples. The following three sample cases 
assume that initial detection is made with northward- 
looking TDWR. Two minutes later, three different events 
are postulated: 



Figure 9: Extent flgures-of-merit: time- 

varying cases 
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time-varying cases 





situation. There is little difference between the three 
sample cases. 

The intensity results (Figure 10) are similar to those 
from the single-microburst time -in variant runs, in that all 
of the estimates are low. As evident from the “TASS 
data" curve, the actual microburst increases in strength in 
the first two minute span and then decreases in the last 
two minutes. Only in scenario 3, in which northbound 
TDWR measurements were combined with eastbound 
ABDR measurements, was the filter able to follow this 
trend. The low estimates are most likely due to 
difficulties matching this complex microburst with the 
single-microburst filter. However, the intensity results 
were somewhat sensitive to the choice of process noise 
strength, which indicates a need for further study of 
microburst time dynamics. 

Discussion 

The simulations demonstrate the potential usefulness 
of this technique, particularly for estimating the size and 
position of the microburst hazard region. Several other 
characteristics of the algorithm were also observed during 
the simulation runs, although it should be noted that the 
use of computational data for a single historical 
microburst event limits the scope of the conclusions that 
can be drawn. Further simulation work is planned, using 
data from actual field measurements. 

The single -microburst algorithm appeared to be 
numerically robust. Errors in initial conditions and 
reasonable variations in choice of filter parameters did not 
produce filter instability in either the time-invariant or the 
time-varying simulations. The two-microburst form, 
however, was numerically sensitive. In several cases the 
filter diverged during the update iteration, and choice of 
parameters such as the initial covariance matrix appeared 
to have a large impact on the convergence properties of 
the filter. In cases where convergence was reached the 
results tended to be dependent on the actual windfield 
shape. When the windfield had two clear downdraft 
centers, the convergence was steady and the results for 
both intensity and extent were good. In cases where there 
was only one region of high downdraft (such as the first 
time-step of TASS model data) then the two sets of 
microburst parameters either became coincident, or one 
microburst became very weak. This mismodeling 
problem was also apparent in the diagonal elements of the 
covariance matrix; the covariance elements corresponding 
to the unnecessary microburst parameters grew very large. 
Possible solutions to this problem include more 
intelligent initialization based on recognized windfield 
features, or running multiple filters of different types in 
parallel. In any case, the improved estimation possible 
from the two-microburst filter must be weighed against 
the associated numerical difficulties. 

Aside from numerical robustness and algorithm 
tuning issues, there are other implementation issues to be 
considered. The computational requirements of the filter 
need to be assessed with respect to available 


computational resources. Computational load can be 
decreased by thinning large data sets, at the expense of 
estimation accuracy. Also, dataiink bandwidth needs to be 
considered. A model parameter list of 9 elements, for 
example, has an associated 81 element covariance matrix 
(of which only 45 are unique). It is likely that the entire 
covariance matrix is not necessary to initialize the next 
update step, and that some elements could be omitted 
without loss of performance. 

Although the algorithm has been presented in the 
context of multi-sensor data fusion, it does not require 
multiple sensors. Benefits would still be gained if it were 
used with a single sensor due to the additional information 
contained in the analytical model (correct fluid dynamics, 
empirical data). Also, the algorithm could be adapted to 
other fluid dynamic phenomena which can be represented 
by simple analytical models. 

Conclusion 

A recursive model-based data fusion algorithm for 
multi-sensor microburst hazard assessment was presented. 
A simple analytical microburst model is used to 
approximate the actual windfield, and a "best" set of model, 
parameters are estimated from measured winds using Ih 
extended Kalman filtering technique. The resulting 
parameter estimate and associated error covariance 
encapsulate the current state of knowledge about the actual 
windfield, and can be used to compute estimates of 
microburst position, extent, and intensity for alert 
generation. Microburst state dynamics and process noise 
parameters for the filter were chosen based on statistical 
data from microburst field studies. 

Simulated measurements for three types of sensors 
were derived from a time- varying computational model of 
a historical microburst event. Two forms of the 
algorithm were then tested, one using a single-microburst 
model and one using a two-microburst model. It was 
found through both time- invariant and time-varying 
simulations that both forms of the algorithm were able to 
estimate the position and extent of the simulated 
microburst well. The two-microburst model produced 
better intensity estimates, but suffered from numerical 
robustness problems. These preliminary results are 
promising, and further work is planned including 
simulations using field measurements and study of 
feasibility issues such as computational requirements. 
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Appendix A: Oseguera-Bo wles-Vicrov 
Microburst Model Equations 


The U, V, and W wind components are functions of 
position (x,y,z) and the model parameters: 
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where the position offsets are given by 
x = x - Xo 
y = y-yo 

and the radial scale factor X is: 
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Ci and C 2 are empirically adjusted constants with the 
following values: 


Ci = -0.15 (A7) 

C 2 = - 3.2175 (A8) 

and a is a shaping parameter which was set to 2.0 for 
the work presented here. 
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